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ABSTRACT 



Context. Supersonic nonthermal motions in molecular clouds are often interpreted as long-lived magnetohydrodynamic (MHD) 
waves. The propagation and amplitude of these waves is affected by local physical characteristics, most importantly the gas den- 
sity and the ionization fraction. 

Aims. We study the propagation, reflection and dissipation of Alfven waves in molecular clouds deriving the behavior of observable 
quantities such as the amplitudes of velocity fluctuations and the rate of energy dissipation. 

Methods. We formulated the problem in terms of Elsasser variables for transverse MHD waves propagating in a one-dimensional 
inhomogeneous medium, including the dissipation due to collisions between ions and neutrals and to a nonlinear turbulent cascade 
treated in a phenomenological way. We considered both steady-state and time-dependent situations and solved the equations of the 
problem numerically with an iterative method and a Lax-Wendroff scheme, respectively. 

Results. Alfven waves incident on overdense regions with density profiles typical of cloud cores embedded in a diffuse gas suffer 
enhanced reflection in the regions of the steepest density gradient, and strong dissipation in the core's interior These effects are es- 
pecially significant when the wavelength is intermediate between the critical wavelength for propagation and the typical scale of the 
density gradient. For larger wave amplitudes and/or steeper input spectra, the effects of the perpendicular turbulent cascade result in 
a stronger energy dissipation in the regions immediately surrounding the dense core. 

Conclusions. The results may help to interpret the sharp decrease of line width observed in the environments of low-mass cloud cores 
in several molecular transitions. 
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1. Introduction 

The supersonic nonthermal motions observed in molecular 
clouds have been attributed to the presence of magnetohydro- 
dynamic (MHD) waves of sub-Alfvenic amplitudes (Arons & 
Max 1975), of which the Alfven waves are the longest-lived 
(Zweibel & Josafatsson 1983). Indeed, molecular clouds and the 
cores embedded in them are threaded by the Galactic magnetic 
field, as revealed by polarization maps of the thermal dust emis- 
sion (see e.g. Tang et al. 2009). 

Molecular cloud cores, the sites of star formation, are char- 
acterized by velocity dispersions with nonthermal motions sub- 
sonic and uniform, in contrast to the surrounding lower density 
gas, where nonthermal motions are supersonic and in general 
increasing with size (Barranco & Goodman 1998; Goodman et 
al. 1998; CaselH et al. 2002). The NH3 observations by Pineda 
et al. (2010) of the B5 region in Perseus revealed that the tran- 
sition in velocity dispersion between the quiescent core and its 
more "turbulent" environment is sharp, taking place in a thin 
layer with a width of less than 0.04 pc. 

Several groups have performed simulations and developed 
(semi)analytical models of the propagation of non-linear and 
linear Alfven waves in interstellar clouds, but the main focus 
of these studies was either on the generation of density fluctu- 
ations and velocity dispersions with properties comparable to 
those of molecular clouds (Heitsch & Burkert 2002; Kudoh & 
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Basu 2003, 2006), or on the role of internally generated Alfven 
waves to support clumps within molecular clouds against their 
own self-gravity, taking into account the ionization structure and 
the ion-neuti-al drift (Martin et al. 1997; Coker et al. 2000). Our 
main objective is to study how density gradients, ionization frac- 
tion and non-linear wave interaction affect the propagation and 
the amplitude of Alfven waves at the small scale of the interface 
between a dense core and its environment. 

For this goal, we considered a simple one-dimensional cloud 
model for which the propagation of Alfven waves is formu- 
lated in terms of Elsasser variables. This formulation allows 
one to treat propagation, reflection, damping, and nonlinear 
interactions in a simpler form compared to the original fluid 
equations (Elsasser 1950; Dobrowolny et al. 1980; Marsch & 
Mangeney 1987). Our approach is similar to that adopted for 
studying the propagation and reflection of Alven waves in radi- 
ally stratified stellar atmospheres and winds (Velli 1993). In this 
context the density gradients have the main effect of modulat- 
ing the amplitude of Alfven waves. This may affect the damping 
that results from the ion-neutral drift and at the same time trigger 
the development of a (perpendicular) turbulent cascade. Here we 
will focus mostly on the first aspect, although we will also use 
a phenomenological model to describe the dissipation of turbu- 
lence (successfully applied to the solar wind, e.g. Verdini et al. 
2010) to see how it affects the results. 

The paper is organized as follows: in sect. |2]we derive the 
wave equation and review the conditions for wave propagation; 
in sect. [3] we formulate the problem in terms of Elssasser vari- 
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ables; in sect. |4]we describe the method of solution; we discuss 
some special cases in sect. |5] and, in sect. |6l we consider wave 
propagation, damping and reflection in a one-dimensional cloud 
model with a specific density and ionization profile; we evaluate 
the energy dissipation due to both ion-neutral drift and nonlinear 
effects; finally in sect.|7]we summarize our conclusions. 



2. Wave equation and dispersion relation 

We considered a one-dimensional density distribution p(x), 
threaded by a magnetic field B - Bx, in which the field is sta- 
tionary and uniform, the latter condition resulting from the re- 
quirement of zero divergence, V ■ B = 0. The total density of the 
gas p is the sum of the densities of neutrals (p„), ions (p,) and 
electrons (pe), and we assumed that the gas is electrically neutral. 
For applications to molecular clouds, we assumed that the field is 
frozen in the ions, which transfer momentum to the neutrals via 
elastic collisions. We neglected the electron's inertia and focused 
on the ambipolar diffusion regime (see e.g. Pinto et al. 2008), ig- 
noring Hall and Ohmic diffusion. In the regime of interest in this 
work («(H2) ~ 10'*-105 cm-\ B ~ 10-100 fiG), Ohmic diffu- 
sion is about 10 orders of magnitudes smaller than ambipolar 
diffusion and can be safely neglected; on the other hand. Hall 
diffusion is about one order of magnitude smaller than ambipo- 
lar diffusion, for standard grain size distribution, and about 4 
orders of magnitudes smaller in the case of large (micron size) 
dust grains (see, e.g., Wardle & Ng 1999; Nakano et al. 2002; 
Pinto et. al. 2008). The consequences of the Hall effect on the 
propagation of Alfven wawes in non-homogeneous media have 
been considered by Campos & Isaeva (1993). 

In our analysis, we also assumed that the magnetic force is 
transferred to the neutrals only by collisions with ions. Actually, 
in the density regime of interest, charged dust grains can con- 
tribute as much as ions in the frictional force, provided that 
a population of small (< fj.m size) grains is present (Nakano 
et al. 2002). If only large (> jum size) grains are present, as 
expected if some amount of dust coagulation has taken place 
(Ormel et al. 2009), the contribution of grains to the frictional 
force becomes much smaller On the other hand, large-size dust 
grains may induce dispersive effects on the propagation of non- 
linear magnetic disturbances, and even promote the formation of 
soliton structures, although on scales at best one order of mag- 
nitude smaller than the typical size of a molecular cloud core 
(Pandey et al. 2008). 

We considered transverse disturbances that are characterized 
by magnetic field b and velocities u„ and u,- (the latter were also 
assumed to be incompressible). We assumed that the sources 
of the waves are located a sufficiently large distance from the 
core's midplane to neglect their density-pressure feedback on 
the cloud. The equations of motion and the induction equation 
read 

Pn-^ YmPiPniUi ' U„) = -p„U„ ■ Vj_U„, (1) 

Ot 



dui B db , ^ „ 1 ,^,2 

Pi-r- - -r-T ~ yinPiPniUn ' u,) = -p,u,- ■ V^u,- - — Vi|b| 
ot 4-7T ox Sn 

+ ^b-V^b, (2) 



where 



{crw)i„ 
Mi + m„ ' 



(4) 



is the coUisional drag coefficient, {a-w)i„ is the momentum trans- 
fer rate coefficient, and is the gradient along the direction 
perpendicular to the large-scale field B. On the right-hand side 
we have grouped the nonlinear terms: note that because of the as- 
sumption of transverse fluctuations they involve only gradients 
in that direction. 

Combining the linearized Eqs. ([U, ^ and (O, we obtain the 
wave equation for u,. 



Ai 



dx^dt 



+ JinPn 



P d^Uj 

Pn dfi 



An 



5V 



= 0, 



(5) 



where p - p„ + p, and Vai, Vau are the Alfven velocities in the 
ions and the neutrals, respectively: 



Vai 



B 



B 



y/^^Pn 



(6) 



Substituting in Eq. (|5]l a solution in the form Ui{x, f) - 
u,o exp[/(A:x- wf)]i we obtain the well- known dispersion relation 
(DR) 



a){a? - {kVAif] + ijinipu" -piikVAif] = 0, 



(7) 



first derived by Piddington (1956) and later by Kulsrud & 
Pearce(1969). 

The amplitudes of the fluctuations in u„ and u, are related to 
the amplitude of the fluctuations of b by 



uo„ = -Va;- 



ijinPi I kVAi 



+ IJinP \ W 



bo 

B ' 



and 



U()/ 



-VAi 



aj[a) + ijinp] - ylfiiPn ( kVAi 
[oj + iyi„p](a> + ijinPn) 



bo 
B ■ 



(8) 



(9) 



For pi/p„ «c 1, there are two separated regimes of wave 
propagation: long waves (for kVAi <K yjpil PnJinPn) and short 
waves (for kVAi » yi„p„). The physical reason is that when 
the timescale for wave damping becomes shorter than the wave 
period, wave oscillation cannot occur. Short waves are high- 
frequency waves propagating at the ion's Alfven speed, damped 
by collisions with neutrals. In the short-wave regime, the ratio of 
the amplitudes of the velocity oscillations reduces to 



lUOnI 

|uo/l 



JinPi 
(X) 



(10) 



i.e. the neutrals are weakly coupled to the ions and their ampli- 
tude of oscillation is very small. Therefore, these waves cannot 
be responsible for the linewidths of neutral species observed in 
molecular clouds. 

In the limit of long waves, the wave equation (|5]l becomes 



5V 



- -r-r - T]AD 



/2[ 



dx^dt 



= 0, 



(11) 



where Va is the Alfven speed in the gas of neutrals and ions (of 
density p - Pn +Pi), and jjad is the ambipolar diffusion resistivity 



(9b du B 
--5^.b.V.u,-u,.V.b, (3) 



(12) 
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In this regime, the DR of Eq. ^ reduces to 



0. 



(13) 



For the propagation and damping of Alfven waves in a molec- 
ular cloud a description in terms of forced waves is appropri- 
ate. Indeed, there are several possible sources of hydromagnetic 
waves in a molecular cloud, such as gravitational collapse of 
smaller parts, large-scale non-radial oscillations of the cloud, ra- 
diation pressure from young stars, expansion of HII regions, stel- 
lar outflows (Arons & Max 1975). In this case, the DR should be 
considered as a quadratic equation for the wavenumber k (with 
real and imaginary part) as a function of the frequency a> (real). 
Defining a complex wave number k - %(k) + i3{k), in the limit 
CO <K V^/rjAD, we derive from the real and imaginary parts of the 
dispersion relation Eq. ( fT3] ) the expressions 



and 



5ik) ^ ± 



2V 



3'7ad- 



(14) 



The condition for wave propagation, |^R(A:)| > |3(^)|, is satisfied 
only if the frequency of forcing is lower than a critical frequency 
LJcr = '^V^/rjAD, or, the wavelength is longer than a critical value 
/Ici - jttiad/Va (Kulsrud & Pearce 1969). In the low-ionization 
regime typical of a molecular cloud, p x p„, the condition for 
wave propagation reduces to <y < 2yi„pj or, in terms of wave- 
length, to 



A> Ar 



o 1/2 

2rmP« p. 



JinPi 



(15) 



We selected the case of depleted clouds where all metal species 
are frozen onto grains and the dominant ions are H"^ and Hj . In 
this case y,-,, = 2.38 x 10"* g"' cm"* s"' (see Pinto & Galli 2008) 
and the mean mass of neutrals and ions are m„ = 2.33mH and 
m, - 3mH, respectively. Inserting numerical values, we obtain 

icr = 3.8 X 10"^^ I ttttt-tt] i — ^'^ (t^-t) ' pc, (16) 

\ 100 jjG I \ 105 cm-^ / \ 10-^ ; F . V 7 

where x, = «, /« is the ionization fraction. In the regime of long 
waves, the amplitudes of the velocity oscillations reduce to 



uo« 



' CD B' 



and 



wbo 



(17) 



Since \k\ - w/Va to lowest order in (?7ad'jJ/V'^) , the two ve- 
locity amplitudes are equal in modulus |uo„| = |uo,| = V/ilbol/B. 
These long-wavelength Alfven waves propagate with speed Va 
in a well-coupled ion-neutral gas, and may be responsible for 
the observed linewidths of ionic and neutral species observed in 
molecular clouds. 



3. Model equations 

The condition for the propagation of long waves is equivalent to 
assume that neutrals and ions are well coupled (|u, - u„|/Va;i ^ 
|b|/B). With this condition, and the approximation p, <K p„, the 
ion velocity can be obtained from Eq. Q as 

'7AD C^b TIkD, _ , ^AD_ ,.,2 ,,o\ 

Upon substitution of the above equation into the induction 
Eq. (O one obtains a quite complicated equation. To make it 
more tractable, still accounting for non-linear interactions, we 



adopted the following additional assumptions: i) fluctuations 
have small amplitudes {bIB ^ u„/Va < 1), thus, when substitut- 
ing Eq. ( fTSl l into Eq. Q we expand the non-linear terms only to 
the second order (cubic terms 0{b^) are neglected); ii) the am- 
bipolar diffusion coefficient varies on a scale much larger than 
all the other quantities do, thus we neglect terms involving its 
derivatives. Under these assumptions, the induction equation is 
of the form 



db du„ 

= B— + b ■ V^u,, 
of ox 



u„ V^b 



(19) 



+J]AD 



d^yd^j'^B ' ^[d^l'lB dx 



We see that nonlinearities also appear in terms that contain the 
ambipolar diffusion coefficient (the last two terms), therefore 
we expect that the turbulent cascade is also affected by ion- 
neutral drift. In view of the simple phenomenology that we 
will use (Eq. ( |26] |) to describe the dissipation of turbulence, 
we additionally simplified the equation by adopting an order- 
ing of the parallel and perpendicular length scales, namely iii) 
Acrlt\\ « b/B « where £\\^ j_ are the typical sizes of 

parallel and perpendicular large scales. We thus obtain the final 
model equations: 



5u„ B db 1 



1 



^^TTT- ^'^^Ibr + ^b ■ V^b -p„u„ ■ V^u„, (20) 
of 4-7T ox isn 4-n 



and 

db du„ d (db\ 

= B— + + b ■ ViU„ - u„ ■ Vib. 

of ox ox \ox 



(21) 



Eq. (l20l i and ( ITTT i can be conveniently expressed in terms of 
the Elsasser variables. 



z+ = u„ + 



y47rp„ 



(22) 



by dividing Eq. ( 1211 1 by ^Anp„, then adding and subtracting 
Eq. (|23 and (HI]): 

^ ± Va$^ = -Z^ ■ V^Z± - —V^P,n ± -^(Zi - Z^) 

of ox p„ 2 



J] AD , . 

±— (z±-Zt) 











5z^\ 


K0- 






, dx 


dx j 



HAD I dhj, dh:f 



2 \ dx^ dx^ 



(23) 



The superscripts ' and " indicate first- and second-order spatial 
derivatives along the x direction and is the magnetic pressure. 
The first two terms on the left-hand side represent the temporal 
and spatial wave propagation, respectively . On the right-hand 
side, the first two term account for non-linear interactions; the 
term proportional to describes the losses related to the pres- 
ence of gradients in the Alfven velocity and wave reflection; the 
terms proportional to 77AD are associated to the dissipation by 
ambipolar diffusion and combined effects of density gradients 
and ion-neutral drift. Because we assumed fluctuations to lie in 
the plane perpendicular to B, the gradient involves only per- 
pendicular scales. Given the high Reynolds numbers found in 
molecular clouds {R„, x 10'^, Myers & Khersonsky 1995), a per- 
pendicular turbulent cascade is expected to develop. These non- 
linear couplings involve only counter-propagating waves. Note 
however that even if only one mode is present, say z+, density 
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gradients or ion-neutral drift ensure the presence of the other 
mode (z-) and hence the triggering of the cascade. 

The hypotheses that we have adopted in deriving 
Eqs. ( I20ll2ni23l l are summarized in the following ordering: 



« — « — « — . 

B h Pi 



(24) 



Note that the first inequality is a stronger constraint than the 
long wavelength condition Eq. ( fTSl l, but in molecular clouds one 
typically finds that Acrli\\ ~ i]ad/VaL < 0.06 (see Eq. (fT6T l) 
and p„/pi > 10^, therefore the above ordering is generally sat- 
isfied even for > 1. The condition b/B « can 
be rewritten as taIt^l < 1 where the Alfven and non-linear 
timescales are given by ta - ^||/Va and t^l - i^lz± (we as- 
sumed bjB «i uIVa ~ z±IVa)- In other words, the equations are 
suited to describe weak turbulence. Our inequality is more strin- 
gent than taItnl ^ Bib given in Galtier et al. (2002), which 
was used to build a cascade model. We will not attempt such a 
detailed description of the cascade for two reasons that are re- 
lated to the inhomogeneity. It will be clear in the following that 
the density gradients couple counter-propagating waves. A first 
consequence is that the phase of oppositely propagating low- 
frequency waves cannot be considered as random anymore, pre- 
venting us from exploiting a phase averaging method at large 
scales to obtain the kinetic equations. Moreover, in the inho- 
mogeneous case typical of the cloud-core systems, the density 
gradients act as semi-reflecting boundaries, which enhance the 
amplitudes, and therefore the cascade rate of large-scale modes. 
These low-frequency waves in which energy is accumulated 
pose a problem with the continuity of the low-frequency tail of 
the spectrum. Our nonlinear model is really only in the spirit 
of von Karman and Howarth (1938), (see also Mac Low 1999), 
where the large eddies control the dissipation rate and the spe- 
cific mechanism for small-scale dissipation is not of central im- 
portance. Large-scale fluctuations control the cascade, and fast, 
small-scale dissipation mechanisms, whether collisional or not, 
act as passive absorbers of cascaded energy. 
Another comment concerns the possibility of nonlinear com- 
pressible interactions and processes such as parametric decay. 
These interactions may play a role, but in a completely tur- 
bulent situations where fluctuations have completely random 
phases, the decay via these resonant processes is greatly reduced 
(Malara & Velli 1996). Consequently, the density fluctuations 
that are there will play a role via processes such as phase mix- 
ing/resonant absorption, and shock heating, but will not alter, ex- 
cept at shocks, the dominant nature of the energy cascade given 
by the dimensional model described below. 



4. Numerical method 

We first sought solutions in form of monochromatic waves for 
the linearized Eqs. ( |23] |. This allowed us to show the role of in- 
terference between waves as produced by reflection and ambipo- 
lar diffusion (see section|5]l. The numerical method is presented 
in the next section 14711 Then we solved the full set of Eqs. ( l23T l 
in their time-dependent form, treating the non-linear interactions 
in a phenomenological way and considering fluctuations that are 
characterized by a steep frequency spectrum. The details of the 
numerical method and the phenomenology adopted are given in 
section r 



4.1. Time-independent sclieme for linear equations 

We neglected non-linear terms and hence the perpendicular 
structure of Alfven waves. We sought solutions for monochro- 
matic waves of the form Z-t(x,t) - VAoy±ix)e'x.p(-i(L)t). For a 
cloud of typical size L and density po, we can define a character- 
istic Alfven speed Vao - B/ -sjA-npo and Alfven time Ia - L/Vao, 
and the nondimensional variables ^ = x/L and p = p/po- In 
nondimensional form, eqs. (l23T l can be written as 



+?7AD^(yi-3';) + ^(3'i'-3';'), 



p" 



(25) 



where W - ojL/Vao and^AD - t1ad/(VaqL). For uniform density 
(p = 1) and without dissipation (^ad - 0), Eq. ( |25] | describes the 
propagation of Alfven waves traveling to the right (+ sign) and 
to the left (- sign). 

Because ambipolar diffusion is a small effect in typical ISM 
conditions (^ad ^ 1), we adopted a perturbational approach. 
We initially solved Eqs. ( l25T l guessing at a functional form for the 
second-order derivatives. In particular, we assumed the exponen- 
tially damped behavior predicted by analytical studies (we have 
verified that the numerical results are independent of the choice 
of the initial guess); then we used the second-order derivatives 
of the solutions obtained by the guess function to integrate 
Eqs. ( |25] | and obtained a new set of second-order derivatives. 
This procedure was iterated until convergence of the solution 
was attained. 

Special care has to be taken when choosing boundary con- 
ditions suitable to describe waves incident on both sides of an 
interstellar cloud because of the interference between the incom- 
ing waves. When only ambipolar diffusion is present, the itera- 
tive procedure reduces the problem to a boundary- value problem 
for the fields z"^. For wave propagation in a density gradient with- 
out dissipation, the boundary-value problem is well defined only 
when one wave vanishes at one boundary. In the more realistic 
case, where both modes are injected at the boundaries, solutions 
depend on the relative phase between the injected z±. The phase 
in fact controls the interference between the incident wave and 
the reflected wave inside the density profile. We will describe 
in detail the appropriate initial conditions in the application to 
interstellar clouds (sect. |6]l, first dealing with density gradients, 
then considering both ambipolar diffusion and density gradients. 

4.2. Time-dependent sclieme for non-linear equations 

To evaluate the effect of turbulent dissipation, which can be trig- 
gered by non-linear wave interactions, we did not model the tur- 
bulent cascade. Instead we used a phenomenological expression 
(Dmitruk et al. 2001) that accounts for the damping of the waves 
and preserves the structure of non-linear terms, allowing interac- 
tions only among counter-propagating waves. We therefore inte- 
grated eqs. ( |23] ) modified with the following substitution in the 
RHS, 



1 

Pn 



|Zt|z± 



(26) 



where {j_ represents the dimension of the largest eddies (the en- 
ergy containing scale in the turbulent field). 

The form of the nonlinear term may be heuristically derived 
in the framework of reduced MHD in which variations along 
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the perpendicular directions are decoupled from those along the 
magnetic field (V = + Vy, with » Vy). In presence of a 
strong magnetic field, the propagation time of the Alfven waves 
Ta is equal to or shorter than, the characteristic time-scale for 
nonlinear interaction t^i over most of the Fourier space (con- 
sistent with the ordering in Eq. (l24ll). The nonlinear cascade 
becomes anisotropic, developing preferentially in the plane per- 
pendicularto the direction of the mean field (Shebalin et al 1983; 
Ougthon et al. 1994; Goldreich & Sridhar 1995). One can there- 
fore Fourier-decompose Eq. ( l23l l only in (wavevector lying 
in plane perpendicular to the large-scale magnetic field). The 
global effect of this perpendicular cascade can be described by 
means of two quantities at the large scales, namely a common 
similarity (correlation) scale and the average velocity difference 
among points belonging to the same eddy. Identifying these two 
quantities with the integral turbulent length (£j_) and the fluctua- 
tions' amplitude of the Elsasser fields, we can construct a char- 
acteristic timescale t^^ = {±/z±, which accounts for nonlinear 
turbulent interactions in Eq. (l23b . This phenomenology assumes 
that all the interactions are resonant, in practice overestimating 
the heating rate for high- frequency fluctuations, which instead 
would have a reduced nonlinear coupling because of their fast 
decorrelation. On the other hand, the phenomenology still ac- 
counts for the reduction of the cascade and associated damping 
due to the imbalance between the Z-t fields. 

In our approximation, parallel scales are not generated by 
the cascade itself, as is expected to happen for strong turbulence 
(Goldreich & Sridhar 1995). Therefore, here the perpendicular 
turbulent cascade represents a channel for energy dissipation that 
is alternative to ambipolar diffusion (assumed to act only on the 
parallel scale). We finally remark that the ordering in Eq. jTM 
is assumed to be valid at large scales but it breaks at smaller 
scales, where the cascade is expected to become strong, eventu- 
ally reaching the dissipative scale. However, for the purpose of 
estimating the turbulent heating rate, the phenomenology is still 
valid and consistent with the ordering since the overall cascade 
rate will be controlled by the weaker cascade at the modeled 
large scales. 

We integrated the evolution Eqs. ( |23] |. with the substitu- 
tion of Ea.(l26b, using a numerical scheme of first order in time 
and second order in space (a centered Lax-Wendroff scheme). 
Passing to time-dependent equations allows us to better model 
the interference (dependence on the phase) of Alfven waves 
(from incoherent sources) that are assumed to enter from both 
sides of the cloud with random phases. Boundary conditions are 
transparent, no reflection occurs there. Waves are injected into 
the cloud according to 



5z+/5f|,=-z, = g(t), 5z-/5f|,=i = fit). 



(27) 



where the forcing /(f), g{t)) has a spectrum peaked at low fre- 
quencies and declining as cj^^. The escaping waves at the bound- 
aries, dz+ldt at X = L and dz-ldt ?A. x - —L, are not specified 
and obtained as part of the solution. 



5. Special cases 

To gain insight into the nature of the solution, we first consid- 
ered the effects of ambipolar diffusion for waves propagating in 
a uniform medium and then the propagation of waves in a non- 
uniform medium without ambipolar diffusion. 



5.1. Damping by ambipolar diffusion, uniform density 
For a uniform density medium, p - I, Eq. dZSl l reduces to 



+ = iWy± + -fiM)(y'l - yl)- 



(28) 



In the limit of small dissipation, setting e = 1, 
Eq. ( l28b has a general solution of the form 



y± 



(29) 



where = (2eW) ' is the damping length and and a-t are con- 
stants. In dimensional units, 



2V^ 



(30) 



Notice that the damping length can also be expressed as id - 

In this case each wave is damped in the direction of prop- 
agation. The solution contains a primary wave traveling in one 
direction (first term on the RHS), whose damping produces a 
secondary wave traveling in the opposite direction (second term 
on the RHS). Notice that when two primary waves are present 
(fl± + 0), the profile of u„ and b will be determined by their 
interference in a non-trivial way. 

5.2. Wave reflection by density gradients 

We then considered a non-uniform medium where density gra- 
dients imply changes in the Alfven speed, but we neglected am- 
bipolar diffusion. Eq. (IZST i becomes 



+ y'^^iWp"^y 



- P ( 
4p 



(31) 



In terms of the scaled amplitude /t = p^^"'y± and the scaled 
wavevector^ = p'^^W, Eq. ( ISTT i becomes 



f'± = ±iq.f± + 



(32) 



The general properties of this equation have been discussed in 
detail by Velli (1993). Here, it is sufficient to notice that for a 
weak density gradient a Taylor expansion gives q{^) x W + qi^, 
which, inserted into Eq. (l32b . gives the approximate solution 



f± ~ (a± + iPa^)e 



(33) 



where y6 = qijAW^ and a-t are constants. Thus, a density gradi- 
ent produces two effects: a variation of the wave amplitude, as 
expressed by the p''''* scaling (Walen 1944) valid in the WKB 
approximation, and the generation of secondary reflected waves 
(second term on RHS). Notice that this approximate solution is 
similar (waves traveling in opposite directions) to Eq. (|29] |. thus, 
again, when two primary waves are present (o-t + 0) the gen- 
eral solution for u„ and b will depend on their interference in a 
non-trivial way. 

In the absence of damping, the wave amplitude variation is 
governed by the conservation of the energy flux S - S + - S - - 
const, where 5 ± = |/-tp/8 are the fluxes carried by the individual 
modes. For a weak gradient, the amplitudes corresponding to the 
solutions (l33l l are 



|/±r ^al+ AjSa+a^ sin(W^) cos(Wf). 



(34) 
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Thus the density gradient generates reflected waves that induce 
spatial oscillations in the wave amplitude, with half the wave- 
length of the incident wave, but of course at the same time con- 
serves the energy flux (S = (a+-fl^)/8). For a given wavelength, 
the stronger the gradients, the stronger the coupling, and the re- 
flected wave tends to have the same amplitude as the incident 
wave. For a weak gradient the coupling is negligible and the so- 
lution approaches the WKB solution. 

6. Application to molecular cloud cores 

We applied our numerical schemes to the study of the propaga- 
tion and damping of Alfven waves in molecular cloud cores. We 
considered a one-dimensional cloud consisting of a low-mass 
core embedded in a uniform density envelope. For the cloud's 
density we assumed a modified isothermal profile. 



p = lUniHlrienw + (ncoie ' "env) Sech (^)], 



(35) 

with fi = 2.36, Henv = 10"* cm"-', Wcore = 10^ cm""', and ^ = x/L 
where L - 0.04 pc. The slab is threaded by a perpendicular, 
straight and uniform magnetic field of strength 70 fiG and we 
considered the low-ionization regime typical of molecular cloud, 
p ~ p„. We assumed a simple ionization law, 



Pi = Cp\ 

with = 1/2 and C = 4.5 x 10"'^ 



cm 



and r,„ = 2.38 x lO'"* cm^ g"' s"' 
ambipolar diffusion resistivity is 

\2 

;7AD = 3 X 10^" ' 



(36) 

-3/2 these values, 

(Pinto & GalU 2008), the 



B 



100 juG 



I llO^cm^V 



(37) 



The effect of wave reflection is illustrated in the following exam- 
ple, which describes the dependence of the velocity fluctuations 
on the normalized wave frequency W. We neglected dissipation 
(assuming jjad = and {± = oo) and assumed that only one 
outward propagating wave (z+) is injected at the right boundary. 
Solutions are obtained integrating backward from ^ - 3 and do 
not depend on the phase of the injected wave. The resulting ve- 
locity amplitudes are shown in Fig. [T] (large panel) for different 
values of the nondimensional frequency W e [0.6,6]. Clearly, 
for increasing frequency the solution approaches the WKB ana- 
lytical solution and the highest frequency case W = 6 is indis- 
tinguishable from it. The amplitude profiles show variations on 
a scale l/k = Va/cj that changes according to the Alfven speed 
profile (shown in the small panel of Fig.[T]i, being smaller in the 
denser region |^| < 1. This means that small scales are naturally 
formed in the core region. 

The inclusion of diffusive effects requires a careful treatment of 
the boundary conditions. The choice of a symmetric density pro- 
file, with vanishing derivative at the cloud's center and at infinity, 
allows one to match the asymptotic solutions with linear combi- 
nations of the solutions obtained by retaining only the terms that 
involve ambipolar diffusion (see sect. 15.11 1. which we require to 
be exponentially damped waves. In practice, we assumed that at 
X = oo (x - -oo) only leftward (rightward) propagating waves 
are present. The integration was carried out only in one half 
of the domain and we looked for solutions symmetric with re- 
spect to the core's midplane. Since the gradients vanish at x = 0, 
the coupling given by reflection is absent and we can impose as 
boundary condition |z+(0)| = |z-(0)| there. We also averaged over 
eight phase differences (A(p = n7r/4 with n = 0, 1, .., 7) between 
z+(0) and z_(0) to mimic the effect of incoherent wave sources at 




Fig. 1. Large panel: normalized amplitudes of velocity fluctua- 
tions as a function of ^ - x/L for waves incident on the right side 
of a core with the density profile given by Eq. ( |35] ) for various 
values of the nondimensional wave frequency W. No dissipation 
is included. The dashed red line shows the WKB solution. Small 
panel: the red line shows the normalized density profile, and the 
black line the normalized Alfven speed . 



infinity. Results are shown in Fig.|2]for the range of normalized 
frequencies 0.6 <W<b, and over a spatial interval correspond- 
ing to roughly three times the core size. The amplitudes of the 
waves were normalized to the values in ^ = 0. The frequency 
interval explored corresponds approximately to the wavelength 
interval L < A < lOL in the core. In the high-frequency regime 
(W = 6), the damping length is comparable to the scale of the 
density gradient {xj - 1.4 at ^ = 0), and the dominant effect 
is the quasi-exponential damping of the waves on a scale ~ x^. 
On the other hand, in the low-frequency regime {W - 0.6) the 
evolution of the wave amplitude is dominated by the the density 
gradient on a scale ^ ~ 1 , because the damping length is much 
longer than the core's size (xd - 140 at ^ = 0). 

Fig. 13 shows the corresponding profiles of velocity fluctua- 
tions for the same range of frequencies as in Fig. |2] Again an 
average procedure was performed to take into account relative 
phases of uncorrected waves moving in opposite directions from 
distant sources. As in the case of the amplitude of the Elsasser 
variables, the velocity amplitudes are mostly sensitive to the den- 
sity gradient at lower frequencies W ~ 0.6 and mostly sensitive 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 




.5 2.0 2.5 3.0 



Fig. 2. Normalized amplitudes of z- (left panel ) and z+ (right 
panel ) as a function of ^ = x/L in a cloud core with the density 
profile given by Eq. i35[ for values of the nondimensional wave 
frequency from W - 0.6 to W = 6. 
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Fig. 3. Normalized amplitude of velocity fluctuations as a func- 
tion of ^ = x/L in a cloud core with the density profile given by 
Eq. (l35T l for values of the nondimensional wave frequency from 
W = 0.6 to W = 6. 



Fig. 4. Ambipolar diffusion heating as a function of the normal- 
ized distance from the core center For comparison, the red line 
shows the cosmic-ray heating rate for a cosmic -ray ionization 
rate 3 x 10"'^ s"' and energy deposition AQ - 20 eV. 



to the ambipolar damping at the highest frequencies {W ~ 6). 
This is easy to understand given the dependence of the ambipolar 
damping length on the inverse square of the frequency (Ea.[30ll. 
Despite the simplicity and idealization of the model, it is clear 
that the sharp gradients of velocity dispersion observed in the en- 
vironments of dense cores (Barranco & Goodman 1998; Caselli 
et al. 2002; Pineda et al. 2010) cannot be produced by a den- 
sity gradient alone, but require the combined eff'ects of a density 
gradient and ambipolar diffusion. In particular, sharp gradients 
in velocity amplitude are only possible when Xd ~ I, which im- 
plies waves with a frequency such as to give (d ^ L. These waves 
have wavelength A ~ {Aa-Ly^^, intermediate between the critical 
wavelength A^i and the typical scale of the density gradient L. 



the values shown in the figure should be considered as upper 
limits to the actual value of the heating rate. For comparison, 
the red curve shows the cosmic-ray heating rate computed as- 
suming a cosmic-ray ionization rate 3 X 10 s ' and a mean 
energy deposited per ionization AQ - 20 eV (Goldsmith 2001). 
The ambipolar diffusion heating is clearly important both in the 
core and in the envelope, and competes with cosmic rays as a 
heating source for the cloud, especially for the highest values of 
the frequency. The behavior of the heating rate deviates signif- 
icantly from the profile obtained assuming b' ~ kb (see, e.g., 
Hennebelle & Inutsuka 2006), because of the modulation of the 
wave amplitude by the density gradient, and the space depen- 
dence of the damping length. 



6.1. Heating rate 

The energy generated by the friction between streaming par- 
ticles is a significant source of heating for molecular clouds 
(Scalo 1977; Lizano & Shu 1987). The energy produced heats 
the bulk of the gas, increasing its pressure and allowing chemical 
reactions to proceed at a faster rate (Flower, Pineau des Forets & 
Hartquist 1985; Pineau des Forets et al. 1986). Therefore, it is 
important to evaluate the heating rate associated to the dissipa- 
tion of the waves considered in the previous sections. 

In a weakly ionized medium, the rate of energy generation 
by ambipolar diffusion (per unit time and unit volume) is 



H = 



|Bx(VxB)p 



(38) 



(Braginskii 1965, see also Pinto et al. 2008). In the special 
case of low-amplitude Alfven waves propagating in a one- 
dimensional medium, Eq. ( l38T l reduces to 



{H) = ^Ibf , 



(39) 



where the brackets indicate the average over time and the prime 
a space derivative. 

The heating rate as function of position for waves propagat- 
ing in the density profile given by Eq. ( l35b is shown in Fig. 31 as 
before for waves in the frequency range 0.6 < < 6. To show 
the results in physical units, we have assumed b - B. Therefore, 



6.2. Energy dissipation by turbulent cascade 

In the previous section we have considered the dissipation of 
linear Alfven waves caused by ion-neutral collisions. Here we 
estimate the effect of turbulent dissipation as given by the phe- 
nomenological expression ( l26b . We decided to switch off am- 
bipolar diffusion to better explore the effect of turbulence, vary- 
ing its strength through the control parameter taIt^l- Setting 
TlAD = we set Acr - and we did not have to worry about 
the first inequality in Eq. ( l24l i. We imposed at both sides of the 
cloud incoming fluctuations characterized by a frequency spec- 
trum of width IjtA and then falling as uT^: energy is at long 
wavelengths and {\\ ^ VaqIa - L. We here discuss four runs in 
which we varied the turbulent strength ta I t^l by changing the 
input amplitude and perpendicular scale as listed in Table[T] The 
strong turbulence cases are runs A and D, the last even though 
it has a smaller input amplitude. The time-averaged, turbulent 
heating rate per unit volume follows directly from Eq. (|26] | and 
is given by 



P« Pnf(o-c)-^, 



where 



k+P + k-P 



(40) 



(41) 
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Fig. 5. Time-averaged normalized energy fluxes entering the 
cloud from both sides, S + /S = -5) and S -/S -(^ - 5), for 
run B. The normalized density profile is shown for comparison 
by the dotted line. 



RUN 


h/B 




1 TnL 


A 


1 


1 


1 


B 


1/VTo 


1 


1/VTo 


C 


1/10 


1 


1/10 


D 


1/10 


1/10 


1 



Table 1. Nondimensional parameters for the runs: input ampli- 
tude (b/B) and ratio of the perpendicular to parallel scales of 
fluctuations The turbulence strength ta/tnl = b{\\/Bij_ 

is also indicated in the last column. 



is the (nomaUzed) cross helicity, E - \z+\^ + is the total 
turbulent energy, and /(o"c) = -y/l - cr^{ VI + cr^ + VI - cr^) 
is a symmetric function that vanishes at |cr^.| - 1 and increases 
monotonically, reaching a maximum /(O) = 1. 
The small amount of energy at high frequency justifies ne- 
glecting ambipolar diffusion, we verified that {H} « Hj for 
flAD ^ < min(b/B) = 0.1. 

Let us examine the case of intermediate turbulence strength, 

ta/tnl = 1/ VTO (run B). In the linear case with t]ad - 0, the 
coupling between z± in proximity of the density gradients causes 
an enhancement of the respective energy fluxes S ±, which have a 
bump toward the cloud's center This ensures the conservation of 
the net energy flux {S +- S -) which is uniform. When nonlinear 
dissipation is taken into account, the profiles of S + and S _ are 
remarkably different, as shown in Fig. |5] where we also plot the 
density profile for comparison (dotted line). Both energy fluxes 
decrease toward the cloud's center, indicating the damping of the 
wave (heating). In contrast to the linear case, there is a sharp de- 
crease on the core size due to the turbulent heating Ht, which is 
stronger in the whole core (being approximately proportional to 
p„). From Fig.|5]one can also deduce the behavior of ctc, which 
is also equal to (5+ - 5'_)/(5'+ -i- 5'_) (the fluxes are normalized 
to their respective incident value, since 5'+(^ = -5) ~ ^-(^ = 5) 
the normalization factors can be considered equal). The normal- 
ized cross helicity vanishes at the core's center and increases 
rapidly in the core's envelope, maintaing a value Icrd > 0.8 at the 
edges of the core. The cross helicity and the total turbulent en- 
ergy affect the average heating rate per unit mass (Hr/pn), which 
is plotted in Fig.|6]for the four runs (run B corresponding to the 



Fig. 6. Profile of the time-averaged heating rate per unit mass (in 
nondimensional units) for different turbulent strength taItnl - 
1, 1/VlO, 1/10, 1, coiTesponding to runs A, B, C, D as labeled 
on each curve {solid, dashed, dot-dashed and long-dashed lines, 
respectively). The normaUzed density profile is shown for com- 
parison as a dotted line. 

dashed line). For weak turbulence (runs B and C) the profile is 
basically flat, varying by 20% at maximum, and thus reflecting 
Ht ^ Pn- Interestingly, it is flat on a region entirely surrounding 
the core and peaks outside it, showing that the wave damping is 
maximal when passing from the cloud to the core. The relatively 
flat profile can be explained in terms of a compensation of the 
terms related to the cross-helicity and total turbulent energy in 
Eq. ( l40b . Indeed, the cascade is more balanced (cTf. ~ 0) at the 
center, but less energetic, while the turbulent energy grows out- 
side the core, but the cascade is less balanced. For stronger tur- 
bulence (e.g. taItnl = 1) these two effects do not compensate 
each other anymore and the minima and maxima in the heating 
rate are more pronounced {Hj follows the scaling with p„ less 
closely). Note that keeping taItnl - 1 but varying the input 
amplitudes bjB - t^lt\\ - 1/10 the profile does not vary, but 
the overall dissipation decreases (compare runs A and D). One 
can verify that the average dissipation per unit mass scales as 
Hjlpn °^ ^/'''IjL^ while the relative excursion at the core edges 
AHt/Ht grows with tnlI^a- 

7. Conclusions 

We have considered the propagation of Alfven waves in a den- 
sity profile typical of a molecular cloud core embedded in a low- 
density medium. The formulation of the problem in terms of 
Elsasser variables was found particularly appropriate to treat the 
effects of wave reflection and dissipation. We have studied with 
time-dependent as well as time-independent numerical codes the 
interaction of incoherent wave trains incident on the opposite 
edges of a one-dimensional model cloud, determining the am- 
plitude variations associated to the reflection and damping of the 
waves. 

We found that for waves with a wavelength longer than the 
critical wavelength, of the order of a few 10"^ pc in a molecular 
cloud core, but shorter or of the order of the size of the core, a 
significant decrease in the velocity amplitude is produced by the 
combined effects of reflection at the core's boundaries and dissi- 
pation due to ambipolar diffusion in the core's interior For these 
waves, the damping length associated to ion-neutral collisions is 
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of the same order as the size of the core. The resulting behavior 
of the amplitude of velocity fluctuations may explain, at least in 
part, the sharp decrease of line widths observed in the environ- 
ments of low-mass cores (see sect.[T]i. Moreover, for these waves 
the energy generated by ion-neutral drift is a significant source 
of heating both in the core and its envelope (see sect. 16. II ). 

We also considered the turbulent dissipation and found that 
the average heating per unit mass scales as Hjlpn {b/ij_y' - 
T^^, thus becoming significant for high enough wave amplitudes 
(0.1 < b/B < 1). This regime is particularly relevant for molec- 
ular clouds and star-forming regions where a statistical analysis 
of dust polarization maps suggests that the ratio of the strength 
of turbulent magnetic field to the large-scale mean field is of the 
order of 0.1-0.5 (Hildebrand et al. 2009; Houde et al. 2009). We 
used a phenomenological term to mimic the effect of a perpen- 
dicular turbulent cascade. The results of our calculations show 
that the turbulent cascade affects the amplitude of fluctuations 
both at the cloud's edge and in the cloud's interior. In particular, 
the energy dissipation (per unit mass) associated to this process 
is peaked in the regions immediately surrounding the dense cen- 
tral core, indicating a stronger damping of the wave amplitude. 
The relative excursion between maxima and minima AHt/H, 
depends on the turbulence strength t^l/ta, so it can be signif- 
icant even for small amplitudes, provided the turbulent correla- 
tion length is of the order of a few times the size of the core. 
In a more realistic situation, the turbulent cascade, although it is 
stronger in the perpendicular direction, also generates small par- 
allel scales. Therefore we expect that turbulence and ambipolar 
diffusion will be at work at the same time, enhancing the sharp 
variation of velocity amplitude at the interface between the cloud 
and the core. In the present work the investigation of the param- 
eter space was limited by the several assumptions (Eq. (l24l) ) we 
made to obtain the model equations. Interesting features appear 
when we are close to the limit of the allowed parameter space. 
Additional work is needed to address these issues with a more 
realistic and detailed modeling of turbulence that includes the 
effects of ambipolar diffusion. 
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